bhla bhla bhla
# data manipulation and visualization
library(tidyverse)
## ── Attaching core tidyverse packages ──────────────────────── tidyverse 2.0.0 ──
## ✔ dplyr 1.1.4 ✔ readr 2.1.5
## ✔ forcats 1.0.0 ✔ stringr 1.5.1
## ✔ ggplot2 3.5.1 ✔ tibble 3.2.1
## ✔ lubridate 1.9.3 ✔ tidyr 1.3.1
## ✔ purrr 1.0.2
## ── Conflicts ────────────────────────────────────────── tidyverse_conflicts() ──
## ✖ dplyr::filter() masks stats::filter()
## ✖ dplyr::lag() masks stats::lag()
## ℹ Use the conflicted package (<http://conflicted.r-lib.org/>) to force all conflicts to become errors
# correlation plot
library(corrplot)
## corrplot 0.94 loaded
# Functional diversity
library(funspace); library(TPD)
## Loading required package: ks
data_specimens <- read_csv("flmnh_data/Museum Data Updated.csv") |>
mutate(species_updated = paste(Genus, Species, sep = " "),
body.mass = as.numeric(Mass))
## Rows: 284 Columns: 27
## ── Column specification ────────────────────────────────────────────────────────
## Delimiter: ","
## chr (15): Meta-archipelago, Archipelago_group, Site, Genus, Species, Subspec...
## dbl (12): Area, DistanceM, Specimen ID, Wing.Length, Secondary.Length, Kipps...
##
## ℹ Use `spec()` to retrieve the full column specification for this data.
## ℹ Specify the column types or set `show_col_types = FALSE` to quiet this message.
## Warning: There was 1 warning in `mutate()`.
## ℹ In argument: `body.mass = as.numeric(Mass)`.
## Caused by warning:
## ! NAs introduced by coercion
species <- unique(data_specimens$species_updated)
We use the raw data (AVONET_Raw_Data sheet in the Excel
file) in AVONET (Tobias et
al. Ecol Letters 2022), to estimate a mean and SD of trait per
species. Unfortunately, raw data of AVONET does not include the
variation of Body mass, they only included the mean value in their
summary (AVONET2_eBird sheet in the Excel file). We will
pull the other sources of functional traits and estimate mean and SD of
this trait afterwards. We also include the taxonomy information in
AVONET (eBird).
avonet <- read_csv("FunctionalTraits/AVONET_Raw_Data.csv") |>
mutate(species_updated = ifelse(eBird.species.group == "Chlorophonia musica flavifrons",
"Chlorophonia flavifrons",
ifelse(eBird.species.group == "Chlorophonia musica musica",
"Chlorophonia musica",
ifelse(eBird.species.group == "Chlorophonia musica sclateri",
"Chlorophonia sclateri",
Species2_eBird)))) |>
dplyr::select(species_updated, Source, Specimen.number,
Beak.Length_Culmen, Beak.Width, Beak.Depth,
Tarsus.Length, `Hand-wing.Index`, Tail.Length) |>
rename(HWI = `Hand-wing.Index`) |>
filter(species_updated %in% species)
## Rows: 90371 Columns: 26
## ── Column specification ────────────────────────────────────────────────────────
## Delimiter: ","
## chr (11): Avibase.ID, Species1_BirdLife, Species2_eBird, eBird.species.group...
## dbl (13): Data.type, Age, Beak.Length_Culmen, Beak.Length_Nares, Beak.Width,...
## lgl (2): Locality, Country
##
## ℹ Use `spec()` to retrieve the full column specification for this data.
## ℹ Specify the column types or set `show_col_types = FALSE` to quiet this message.
This raw data does not have Mass values. We can extracted from the
summarized AVONET2_eBird dataset.
avonet2_eBird <- read_csv("FunctionalTraits/AVONET2_eBird.csv") |>
dplyr::select(Species2,
Mass, `Hand-Wing.Index`, Tarsus.Length, Tail.Length,
Beak.Length_Culmen, Beak.Width, Beak.Depth) |>
rename(body.mass = Mass) |>
rename(HWI = `Hand-Wing.Index`)
## Rows: 10661 Columns: 31
## ── Column specification ────────────────────────────────────────────────────────
## Delimiter: ","
## chr (13): Species2, Family2, Order2, Avibase.ID2, Mass.Source, Mass.Refs.Oth...
## dbl (18): Total.individuals, Female, Male, Unknown, Complete.measures, Beak....
##
## ℹ Use `spec()` to retrieve the full column specification for this data.
## ℹ Specify the column types or set `show_col_types = FALSE` to quiet this message.
avonet_raw <- avonet |>
left_join(avonet2_eBird |> dplyr::select(Species2, body.mass),
join_by("species_updated" == "Species2"))
Finally, we added the archipelago group (manually) to these 383 specimens
avonet_islands <- read_csv("flmnh_data/AVONET_b_m_specimens_to_review.csv") |>
left_join(avonet_raw)
## New names:
## Rows: 383 Columns: 29
## ── Column specification
## ──────────────────────────────────────────────────────── Delimiter: "," chr
## (13): Avibase.ID, Species1_BirdLife, Species2_eBird, eBird.species.group... dbl
## (14): ...1, Data.type, Age, Beak.Length_Culmen, Beak.Length_Nares, Beak.... lgl
## (2): Locality, Country
## ℹ Use `spec()` to retrieve the full column specification for this data. ℹ
## Specify the column types or set `show_col_types = FALSE` to quiet this message.
## Joining with `by = join_by(Source, Specimen.number, Beak.Length_Culmen,
## Beak.Width, Beak.Depth, Tarsus.Length, Tail.Length, species_updated)`
## • `` -> `...1`
The NA data will be replaced by the mean reported in AVONET
names(data_specimens)
## [1] "Meta-archipelago" "Archipelago_group" "Site"
## [4] "Area" "DistanceM" "Genus"
## [7] "Species" "Subspecies" "Specimen ID"
## [10] "Date" "Data" "Country"
## [13] "Locality" "Elevation" "Other Info"
## [16] "Sex" "Mass" "Wing.Length"
## [19] "Secondary.Length" "Kipps_D" "HWI"
## [22] "Tarsus.Length" "Tail.Length" "Beak.Width"
## [25] "Beak.Depth" "Beak.Length_Culmen" "Measurer"
## [28] "species_updated" "body.mass"
names(avonet_islands)
## [1] "...1" "Avibase.ID" "Species1_BirdLife"
## [4] "Species2_eBird" "eBird.species.group" "Species3_BirdTree"
## [7] "Data.type" "Source" "Specimen.number"
## [10] "Sex" "Age" "Locality"
## [13] "Country_WRI" "Country" "Beak.Length_Culmen"
## [16] "Beak.Length_Nares" "Beak.Width" "Beak.Depth"
## [19] "Tarsus.Length" "Wing.Length" "Kipps.Distance"
## [22] "Secondary1" "Hand-wing.Index" "Tail.Length"
## [25] "Measurer" "Protocol" "Publication"
## [28] "species_updated" "Archipelago_group" "HWI"
## [31] "body.mass"
traits_specimens <- data_specimens |>
dplyr::select(Archipelago_group, species_updated,
body.mass, HWI, Tarsus.Length, Tail.Length,
Beak.Length_Culmen, Beak.Width, Beak.Depth)
traits_avonet <- avonet_islands |>
dplyr::select(Archipelago_group, species_updated,
body.mass, HWI, Tarsus.Length, Tail.Length,
Beak.Length_Culmen, Beak.Width, Beak.Depth)
data_traits <- traits_specimens |>
rbind(traits_avonet) |>
group_by(species_updated) |>
mutate(body.mass = ifelse(is.na(body.mass), avonet2_eBird$body.mass,
body.mass),
HWI = ifelse(is.na(HWI), avonet2_eBird$HWI,
HWI),
Tarsus.Length = ifelse(is.na(Tarsus.Length), avonet2_eBird$Tarsus.Length,
Tarsus.Length),
Tail.Length = ifelse(is.na(Tail.Length), avonet2_eBird$Tail.Length,
Tail.Length),
Beak.Length_Culmen = ifelse(is.na(Beak.Length_Culmen), avonet2_eBird$Beak.Length_Culmen,
Beak.Length_Culmen),
Beak.Width = ifelse(is.na(Beak.Width), avonet2_eBird$Beak.Width,
Beak.Width),
Beak.Depth = ifelse(is.na(Beak.Depth), avonet2_eBird$Beak.Depth,
Beak.Depth))
Correct some names and add Id for each data, grouping by species
colnames(data_traits) <- c("Archipelago", "Species",
"mass", "h.w.i", "tars.l", "tail.l",
"beak.l", "beak.w", "beak.d")
data_traits <- data_traits |>
group_by(Species, Archipelago) |>
mutate(Sequence = row_number()) |>
ungroup()
# a function for extract strings
id_specimen <- function(archipelago, species) {
archipelago_part <- str_extract_all(archipelago,
"\\b\\w{2}") |>
unlist() |>
paste0(collapse = "")
species_part <- str_extract_all(species,
"\\b\\w{3}") |>
unlist() |>
paste0(collapse = "")
paste0(species_part, "_", archipelago_part)
}
data_traits <- data_traits |>
mutate(Group = mapply(id_specimen, Archipelago, Species),
ID = paste0(Group, "_", Sequence))
eup_traits <- data_traits |>
filter(Species != "Coereba flaveola")
eup_traits_scl <- as.data.frame(scale(eup_traits[c(3:9)]))
corrplot(round(cor(eup_traits_scl),2), type="upper", order="hclust",
tl.col="black", tl.srt=45)
rownames(eup_traits_scl) <- eup_traits$ID
# Calculating the dimensionality of a functional space based on PCA
funspaceDim(eup_traits_scl)
##
## Using eigendecomposition of correlation matrix.
##
## Retain the first 3 components
## [1] 3
Run PCA of the scaled traits and building the functional trait spaces
pca.traits.eup <- princomp(eup_traits_scl, cor = T)
plot(pca.traits.eup)
And we can do a functional trait space analysis grouping by species and archipelago
# trait space by groups
trait.space.pca.eu.global <- funspace(x = pca.traits.eup,
n_divisions = 200)
# trait space by Archipelago
trait.space.pca.eu.arch <- funspace(x = pca.traits.eup,
group.vec = eup_traits$Archipelago,
n_divisions = 200)
# trait space by group of species by archipelago
trait.space.pca.eu.group<- funspace(x = pca.traits.eup,
group.vec = eup_traits$Group,
n_divisions = 200)
We can explore the results
summary(trait.space.pca.eu.global)
##
## Functional space based on a PCA with 7 dimensions
## Dimensions 1 and 2 are considered in analyses
##
## Loadings:
## Comp.1 Comp.2 Comp.3 Comp.4 Comp.5 Comp.6 Comp.7
## mass 0.300 0.522 0.226 0.745 0.168 0.046 0.001
## h.w.i -0.636 -0.229 -0.284 0.188 0.115 0.643 -0.018
## tars.l -0.360 0.566 -0.488 0.041 -0.557 -0.018 0.014
## tail.l -0.377 0.657 -0.168 -0.386 0.499 -0.023 -0.015
## beak.l -0.328 -0.335 -0.661 0.335 0.225 -0.424 -0.020
## beak.w -0.859 -0.038 0.387 0.078 -0.013 -0.136 0.294
## beak.d -0.825 -0.007 0.451 0.070 -0.089 -0.151 -0.284
##
## Percentage of variance explained for each trait:
## Comp.1 Comp.2 Overall_explained
## mass 9.03 27.27 36.30
## h.w.i 40.47 5.24 45.71
## tars.l 12.93 31.98 44.92
## tail.l 14.18 43.15 57.33
## beak.l 10.77 11.20 21.97
## beak.w 73.72 0.14 73.86
## beak.d 68.08 0.00 68.08
##
## ------------------------------------------------------------
##
## Functional diversity indicators:
##
## ---> For the global set of species:
##
## Functional richness (99.9% probability threshold) = 47.26
##
## Functional divergence = 0.87
plot(x = trait.space.pca.eu.global, type = "global",
quant.plot = T, arrows = T, arrows.length = 0.75)
summary(trait.space.pca.eu.arch)
##
## Functional space based on a PCA with 7 dimensions
## Dimensions 1 and 2 are considered in analyses
##
## Loadings:
## Comp.1 Comp.2 Comp.3 Comp.4 Comp.5 Comp.6 Comp.7
## mass 0.300 0.522 0.226 0.745 0.168 0.046 0.001
## h.w.i -0.636 -0.229 -0.284 0.188 0.115 0.643 -0.018
## tars.l -0.360 0.566 -0.488 0.041 -0.557 -0.018 0.014
## tail.l -0.377 0.657 -0.168 -0.386 0.499 -0.023 -0.015
## beak.l -0.328 -0.335 -0.661 0.335 0.225 -0.424 -0.020
## beak.w -0.859 -0.038 0.387 0.078 -0.013 -0.136 0.294
## beak.d -0.825 -0.007 0.451 0.070 -0.089 -0.151 -0.284
##
## Percentage of variance explained for each trait:
## Comp.1 Comp.2 Overall_explained
## mass 9.03 27.27 36.30
## h.w.i 40.47 5.24 45.71
## tars.l 12.93 31.98 44.92
## tail.l 14.18 43.15 57.33
## beak.l 10.77 11.20 21.97
## beak.w 73.72 0.14 73.86
## beak.d 68.08 0.00 68.08
##
## ------------------------------------------------------------
##
## Functional diversity indicators:
##
## ---> For the global set of species:
##
## Functional richness (99.9% probability threshold) = 47.26
##
## Functional divergence = 0.87
##
## ---> For each group:
##
## Functional richness (99.9% probability threshold; Continental island) = 5.69
## Functional richness (99.9% probability threshold; Greater Antilles) = 10.35
## Functional richness (99.9% probability threshold; Lesser Antilles - Kalinago) = 3.48
## Functional richness (99.9% probability threshold; Mainland) = 44.03
##
## Functional divergence (Continental island) = 0.76
## Functional divergence (Greater Antilles) = 0.83
## Functional divergence (Lesser Antilles - Kalinago) = 0.39
## Functional divergence (Mainland) = 0.84
plot(x = trait.space.pca.eu.arch, type = "groups",
quant.plot = T,globalContour = T, pnt = T, pnt.cex = 0.1)
summary(trait.space.pca.eu.group)
##
## Functional space based on a PCA with 7 dimensions
## Dimensions 1 and 2 are considered in analyses
##
## Loadings:
## Comp.1 Comp.2 Comp.3 Comp.4 Comp.5 Comp.6 Comp.7
## mass 0.300 0.522 0.226 0.745 0.168 0.046 0.001
## h.w.i -0.636 -0.229 -0.284 0.188 0.115 0.643 -0.018
## tars.l -0.360 0.566 -0.488 0.041 -0.557 -0.018 0.014
## tail.l -0.377 0.657 -0.168 -0.386 0.499 -0.023 -0.015
## beak.l -0.328 -0.335 -0.661 0.335 0.225 -0.424 -0.020
## beak.w -0.859 -0.038 0.387 0.078 -0.013 -0.136 0.294
## beak.d -0.825 -0.007 0.451 0.070 -0.089 -0.151 -0.284
##
## Percentage of variance explained for each trait:
## Comp.1 Comp.2 Overall_explained
## mass 9.03 27.27 36.30
## h.w.i 40.47 5.24 45.71
## tars.l 12.93 31.98 44.92
## tail.l 14.18 43.15 57.33
## beak.l 10.77 11.20 21.97
## beak.w 73.72 0.14 73.86
## beak.d 68.08 0.00 68.08
##
## ------------------------------------------------------------
##
## Functional diversity indicators:
##
## ---> For the global set of species:
##
## Functional richness (99.9% probability threshold) = 47.26
##
## Functional divergence = 0.87
##
## ---> For each group:
##
## Functional richness (99.9% probability threshold; Chlcya_Ma) = 9.65
## Functional richness (99.9% probability threshold; Chlele_Ma) = 18.18
## Functional richness (99.9% probability threshold; Chlfla_LeAnKa) = 3.48
## Functional richness (99.9% probability threshold; Chlmus_GrAn) = 6.11
## Functional richness (99.9% probability threshold; Chlocc_Ma) = 3.74
## Functional richness (99.9% probability threshold; Chlscl_GrAn) = 8.5
## Functional richness (99.9% probability threshold; Eupaff_Ma) = 7.93
## Functional richness (99.9% probability threshold; Eupchl_Ma) = 7.95
## Functional richness (99.9% probability threshold; Eupful_Ma) = 3.13
## Functional richness (99.9% probability threshold; Eupgou_Ma) = 11.34
## Functional richness (99.9% probability threshold; Euphir_Ma) = 3.27
## Functional richness (99.9% probability threshold; Eupjam_GrAn) = 3.98
## Functional richness (99.9% probability threshold; Euplan_Ma) = 30.32
## Functional richness (99.9% probability threshold; Eupmin_Ma) = 6.76
## Functional richness (99.9% probability threshold; Euppec_Ma) = 5.27
## Functional richness (99.9% probability threshold; Eupruf_Ma) = 5.66
## Functional richness (99.9% probability threshold; Eupsat_Ma) = 5.77
## Functional richness (99.9% probability threshold; Euptri_Cois) = 2.72
## Functional richness (99.9% probability threshold; Euptri_Ma) = 3.18
## Functional richness (99.9% probability threshold; Eupvio_Cois) = 5.12
## Functional richness (99.9% probability threshold; Eupvio_Ma) = 13.01
##
## Functional divergence (Chlcya_Ma) = 0.66
## Functional divergence (Chlele_Ma) = 0.77
## Functional divergence (Chlfla_LeAnKa) = 0.39
## Functional divergence (Chlmus_GrAn) = 0.58
## Functional divergence (Chlocc_Ma) = 0.37
## Functional divergence (Chlscl_GrAn) = 0.86
## Functional divergence (Eupaff_Ma) = 0.6
## Functional divergence (Eupchl_Ma) = 0.56
## Functional divergence (Eupful_Ma) = 0.35
## Functional divergence (Eupgou_Ma) = 0.68
## Functional divergence (Euphir_Ma) = 0.39
## Functional divergence (Eupjam_GrAn) = 0.4
## Functional divergence (Euplan_Ma) = 0.78
## Functional divergence (Eupmin_Ma) = 0.6
## Functional divergence (Euppec_Ma) = 0.51
## Functional divergence (Eupruf_Ma) = 0.6
## Functional divergence (Eupsat_Ma) = 0.81
## Functional divergence (Euptri_Cois) = 0.35
## Functional divergence (Euptri_Ma) = 0.38
## Functional divergence (Eupvio_Cois) = 0.82
## Functional divergence (Eupvio_Ma) = 0.83
plot(x = trait.space.pca.eu.group, type = "groups",
quant.plot = T, globalContour = T, pnt = T, pnt.cex = 0.1)
coe_traits <- data_traits |>
filter(Species == "Coereba flaveola")
coe_traits_scl <- as.data.frame(scale(coe_traits[c(3:9)]))
corrplot(round(cor(coe_traits_scl),2), type="upper", order="hclust",
tl.col="black", tl.srt=45)
rownames(coe_traits_scl) <- coe_traits$ID
# Calculating the dimensionality of a functional space based on PCA
funspaceDim(coe_traits_scl)
##
## Using eigendecomposition of correlation matrix.
##
## Retain the first 1 components
## [1] 1
Run PCA of the scaled traits and building the functional trait spaces
pca.traits.coe <- princomp(coe_traits_scl, cor = T)
plot(pca.traits.coe)
And we can do a functional trait space analysis grouping by species and archipelago
# trait space by groups
trait.space.pca.coe.global <- funspace(x = pca.traits.coe,
n_divisions = 200)
# trait space by Archipelago
trait.space.pca.coe.arch <- funspace(x = pca.traits.coe,
group.vec = coe_traits$Archipelago,
n_divisions = 200)
We can explore the results
summary(trait.space.pca.coe.global)
##
## Functional space based on a PCA with 7 dimensions
## Dimensions 1 and 2 are considered in analyses
##
## Loadings:
## Comp.1 Comp.2 Comp.3 Comp.4 Comp.5 Comp.6 Comp.7
## mass 0.237 0.810 0.010 0.451 0.278 0.042 0.072
## h.w.i -0.649 -0.304 -0.015 0.192 0.446 0.500 0.008
## tars.l -0.508 0.307 -0.533 -0.053 -0.518 0.303 0.033
## tail.l -0.308 0.410 0.449 -0.708 0.103 0.151 0.020
## beak.l -0.500 0.091 -0.615 -0.274 0.412 -0.344 -0.028
## beak.w -0.842 -0.068 0.259 0.182 -0.116 -0.251 0.332
## beak.d -0.825 0.145 0.283 0.256 -0.123 -0.160 -0.335
##
## Percentage of variance explained for each trait:
## Comp.1 Comp.2 Overall_explained
## mass 5.62 65.61 71.23
## h.w.i 42.18 9.26 51.44
## tars.l 25.77 9.44 35.21
## tail.l 9.52 16.81 26.33
## beak.l 24.95 0.82 25.78
## beak.w 70.88 0.47 71.35
## beak.d 68.04 2.10 70.14
##
## ------------------------------------------------------------
##
## Functional diversity indicators:
##
## ---> For the global set of species:
##
## Functional richness (99.9% probability threshold) = 39.85
##
## Functional divergence = 0.86
plot(x = trait.space.pca.coe.global, type = "global",
quant.plot = T, arrows = T, arrows.length = 0.75)
summary(trait.space.pca.coe.arch)
##
## Functional space based on a PCA with 7 dimensions
## Dimensions 1 and 2 are considered in analyses
##
## Loadings:
## Comp.1 Comp.2 Comp.3 Comp.4 Comp.5 Comp.6 Comp.7
## mass 0.237 0.810 0.010 0.451 0.278 0.042 0.072
## h.w.i -0.649 -0.304 -0.015 0.192 0.446 0.500 0.008
## tars.l -0.508 0.307 -0.533 -0.053 -0.518 0.303 0.033
## tail.l -0.308 0.410 0.449 -0.708 0.103 0.151 0.020
## beak.l -0.500 0.091 -0.615 -0.274 0.412 -0.344 -0.028
## beak.w -0.842 -0.068 0.259 0.182 -0.116 -0.251 0.332
## beak.d -0.825 0.145 0.283 0.256 -0.123 -0.160 -0.335
##
## Percentage of variance explained for each trait:
## Comp.1 Comp.2 Overall_explained
## mass 5.62 65.61 71.23
## h.w.i 42.18 9.26 51.44
## tars.l 25.77 9.44 35.21
## tail.l 9.52 16.81 26.33
## beak.l 24.95 0.82 25.78
## beak.w 70.88 0.47 71.35
## beak.d 68.04 2.10 70.14
##
## ------------------------------------------------------------
##
## Functional diversity indicators:
##
## ---> For the global set of species:
##
## Functional richness (99.9% probability threshold) = 39.85
##
## Functional divergence = 0.86
##
## ---> For each group:
##
## Functional richness (99.9% probability threshold; Continental island) = 4.55
## Functional richness (99.9% probability threshold; Greater Antilles) = 16.15
## Functional richness (99.9% probability threshold; Lesser Antilles - Kalinago) = 9.24
## Functional richness (99.9% probability threshold; Lucayan) = 11.65
## Functional richness (99.9% probability threshold; Mainland) = 29.17
##
## Functional divergence (Continental island) = 0.34
## Functional divergence (Greater Antilles) = 0.79
## Functional divergence (Lesser Antilles - Kalinago) = 0.55
## Functional divergence (Lucayan) = 0.71
## Functional divergence (Mainland) = 0.82
plot(x = trait.space.pca.coe.arch, type = "groups",
quant.plot = T, globalContour = T, pnt = T, pnt.cex = 0.1)